#ifndef I_LIB_IP_I_BASIC_H
#define I_LIB_IP_I_BASIC_H

#include "../core/i_blas.h"
#include "../core/i_struct.h"

namespace idl
{
    /*This class assumes that the Neighborhood is symmetric, i.e., n is odd NOT even*/
    template<int n, typename T>
    class Neighborhood
    {
        public:
            static const int DIAMETER = n;
            static const int RADIUS = (n >> 1);
            typedef T DATATYPE;
            Neighborhood()
            {
                int i = 0, dr, dc;
                _region_size = n*n;
                //if (n)
                {
                    for (dr = -RADIUS; dr <= RADIUS; dr++)
                    for (dc = -RADIUS; dc <= RADIUS; dc++)
                    {
                        _offsets[i].first = dr;
                        _offsets[i].second = dc;
                        i++;
                    }
                    /*default, uniform weights*/
                    i_fill(_weights, _region_size, i_rec((T)_region_size)); 
                }
            }
            explicit Neighborhood(T sigma)
            {
                int i = 0, dr, dc;
                _region_size = n*n;
                //if (n)
                {
                    for (dr = -RADIUS; dr <= RADIUS; dr++)
                    for (dc = -RADIUS; dc <= RADIUS; dc++)
                    {
                        _offsets[i].first = dr;
                        _offsets[i].second = dc;
                        i++;
                    }
                    if (sigma <= (T)0.0)
                    {
                        /*invalid sigma, uniform weights*/
                        i_fill(_weights, _region_size, i_rec((T)_region_size)); 
                    }
                    else
                    {
                        i_gaussian_2d(_weights, n, sigma);
                    }
                }
            };
        virtual ~Neighborhood(){};
        const T* get_weights() const { return _weights; }
        const Pair<int, int>*  get_offsets() const { return _offsets; }
        int get_size() const { return _region_size; }
    protected:
        int _region_size;
        T _weights[n*n];
        Pair<int, int> _offsets[n*n];/*row offsets*/
    };

    /*This class assumes that the Neighborhood is symmetric, i.e., n is odd NOT even*/
    template<int n, typename T>
    class NeighborhoodCircular : public Neighborhood<n, T>
    {
        public:
        typedef Neighborhood<n, T> Base;
    
        /*for Linux compiler*/
        using Neighborhood<n, T>::RADIUS;
        using Neighborhood<n, T>::_region_size;
        using Neighborhood<n, T>::_weights;
        using Neighborhood<n, T>::_offsets;

        NeighborhoodCircular()
        {
            int i = 0, dr, dc, drdr;
            static const int RR = i_sqr(RADIUS);
            _region_size = 0;
            //if (n)
            {
                for (dr = -RADIUS; dr <= RADIUS; dr++)
                {
                    drdr = i_sqr(dr);
                    for (dc = -RADIUS; dc <= RADIUS; dc++)
                    {
                        if (drdr + i_sqr(dc) <= RR)
                        {
                            _offsets[i].first = dr;
                            _offsets[i].second = dc;
                            i++;
                        }
                    }
                }
                _region_size = i;
                /*default, uniform weights*/
                i_fill(_weights, _region_size, i_rec((T)_region_size)); 
            }
        }
        explicit NeighborhoodCircular(T sigma)
        {
            int i = 0, dr, dc; 
            static const T RR = (T)i_sqr(RADIUS);
            _region_size = 0;
            //const T cen = ((T)(n - 1)) / 2;
            //normalization factor = 1/(2sigma*sigma)
            const T nf = i_div((T)0.5, (sigma*sigma)); 
            T drdr, dcdc, v, ksum = (T)0.0;
            //if (n)
            {
                for (dr = -RADIUS; dr <= RADIUS; dr++)
                {
                    drdr = (T)i_sqr(dr);
                    for (dc = -RADIUS; dc <= RADIUS; dc++)
                    {
                        dcdc = (T)i_sqr(dc);
                        if (drdr + dcdc <= RR)
                        {
                            _offsets[i].first = dr;
                            _offsets[i].second = dc;
                            _weights[i] = v = i_exp(-(drdr + dcdc)*nf);
                            ksum += v;
                            i++;
                        }
                    }
                }
                _region_size = i;
                if (sigma <= (T)0.0)
                {
                    /*invalid sigma, uniform weights*/
                    i_fill(_weights, _region_size, i_rec((T)_region_size)); 
                }
                else
                {
                    /*normalize the kernel*/
                    i_scale(_weights, _region_size, i_rec(ksum));
                }
            }
        };
        virtual ~NeighborhoodCircular(){};
    };

  /*Nearest neighbor sampling*/
    template<typename T>
    inline T i_grid_sampling_nn(const int c[2], const T *const *data, int width, int height)
    {
        T val = (T)0.0;
        int x = c[0];
        int y = c[1];
        if (i_within_2d(x, y, 0, 0, width, height))
        {
            val = data[y][x];
        }
        return val;
    }

    template<typename T>
    inline T i_grid_sampling_nn(const float c[2], const T *const *data, int width, int height)
    {
        int x = i_round(c[0]);
        int y = i_round(c[1]);
        T val = (T)0.0;
        if (i_within_2d(x, y, 0, 0, width, height))
        {
            val = data[y][x];
        }
        return val;
    }

    template<typename T>
    inline T i_grid_sampling_nn(const double c[2], const T *const *data, int width, int height)
    {
        int x = i_round(c[0]);
        int y = i_round(c[1]);
        T val = (T)0.0;
        if (i_within_2d(x, y, 0, 0, width, height))
        {
            val = data[y][x];
        }
        return val;
    }

    /*Bilinear interpolation sampling*/
    template<typename T>
    inline T i_grid_sampling_bilinear(const int c[2], const T *const *data, int width, int height)
    {
        return i_grid_sampling_nn<T>(c, data, width, height);
    }

    template<typename T>
    inline T i_grid_sampling_bilinear(const float c[2], const T *const *data, int width, int height)
    {
        float u = c[0];
        float v = c[1];
        int iu = (int)u;
        int iv = (int)v;
        int iup, ivp;
        float ks, kb, tmpv1, tmpv2;
        T val = (T)0.0;
        if (i_within_2d(iu, iv, 0, 0, width, height))
        {
            iup = i_min(iu + 1, width - 1);
            ivp = i_min(iv + 1, height - 1);
            ks = (float)(iup) - u;
            kb = 1.0f - ks;
            tmpv1 = (ks * (float)data[iv][iu]) + (kb * (float)data[iv][iup]);
            tmpv2 = (ks * (float)data[ivp][iu]) + (kb * (float)data[ivp][iup]);
            ks = (float)(ivp) - v;
            kb = 1.0f - ks;
            val = (T)(ks * tmpv1 + kb * tmpv2);
        }
        return val;
    }

    template<typename T>
    inline T i_grid_sampling_bilinear(const double c[2], 
        const T *const *data, 
        int width, int height)
    {
        double u = c[0];
        double v = c[1];
        int iu = (int)u;
        int iv = (int)v;
        int iup, ivp;
        double ks, kb, tmpv1, tmpv2;
        T val = (T)0.0;
        if (i_within_2d(iu, iv, 0, 0, width, height))
        {
            iup = i_min(iu + 1, width - 1);
            ivp = i_min(iv + 1, height - 1);
            ks = (double)iup-u;
            kb = 1.0 - ks;
            tmpv1 = (ks * (double)data[iv][iu]) + (kb * (double)data[iv][iup]);
            tmpv2 = (ks * (double)data[ivp][iu]) + (kb * (double)data[ivp][iup]);
            ks = (double)(ivp)-v;
            kb = 1.0 - ks;
            val = (T)(ks * tmpv1 + kb * tmpv2);
        }
        return val;
    }

    inline void i_lum8_to_rgb8(const unsigned char* lum8, unsigned char* rgb8, int size)
    {
        int i, i3;
        for (i = 0, i3 = 0; i < size; ++i, i3 += 3)
        {
            i_fill3(rgb8 + i3, lum8[i]);
        }
    }

    inline void i_lum8_to_lum32f(const unsigned char* lum8, float* lum32f, int size)
    {
        int i;
        for (i = 0; i < size; ++i)
        {
            lum32f[i] = (float)(lum8[i]);
        }
    }

    inline void i_lum8_to_lum32f(const unsigned char* lum8, 
        float* lum32f, 
        int size, 
        float scale_factor)
    {
        int i;
        for (i = 0; i < size; ++i)
        {
            lum32f[i] = (float)(lum8[i]) * scale_factor;
        }
    }

    inline void i_rgb8_to_lum8(const unsigned char* rgb, 
        unsigned char* lum, 
        int size)
    {
        int i, j;
        for (i = 0, j = 0; i < size; ++i, j += 3)
        {
            lum[i] = (unsigned char)i_round(((float)rgb[j] * 0.299f) 
                + ((float)rgb[j + 1] * 0.587f) 
                + ((float)rgb[j + 2] * 0.114f));
        }
    }

    inline void i_bgr8_to_lum8(const unsigned char* rgb, unsigned char* lum, int size)
    {
        int i, j;
        for (i = 0, j = 0; i < size; ++i, j += 3)
        {
            lum[i] = (unsigned char)i_round(((float)rgb[j] * 0.114f) 
                + ((float)rgb[j + 1] * 0.587f) 
                + ((float)rgb[j + 2] * 0.299f));
        }
    }

    inline void i_rgb8_to_lum32f(const unsigned char* rgb, float* lum, int size)
    {
        int i, j;
        for (i = 0, j = 0; i < size; ++i, j += 3)
        {
            lum[i] = ((float)rgb[j] * 0.299f) 
                + ((float)rgb[j+1] * 0.587f) 
                + ((float)rgb[j+2] * 0.114f);
        }
    }

    inline void i_bgr8_to_lum32f(const unsigned char* rgb, float* lum, int size)
    {
        int i, j;
        for (i = 0, j = 0; i < size; ++i, j += 3)
        {
            lum[i] = ((float)rgb[j] * 0.114f) 
                + ((float)rgb[j + 1] * 0.587f) 
                + ((float)rgb[j + 2] * 0.299f);
        }
    }
} /*namespace idl*/

#endif